pacman::p_load(ggrepel, patchwork,
ggthemes, hrbrthemes,
ggdist, ggridges,
colorspace, gridExtra,
tidyverse, tmap, sf, tmaptools, dplyr, tibble, dplyr) Big Mac Index Choropleth and Time Series
Take Home Exercise 04
1 Load Packages and Data
Importing data
big_mac <- read_csv("data/countries_with_complete_data.csv")2 ChoroPleth
data("World")wmpb_World <- World %>%
select(-HPI, -inequality, -footprint, -life_exp, -well_being, -gdp_cap_est)
wmpb_WorldSimple feature collection with 177 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
Geodetic CRS: WGS 84
First 10 features:
iso_a3 name sovereignt continent
1 AFG Afghanistan Afghanistan Asia
2 AGO Angola Angola Africa
3 ALB Albania Albania Europe
4 ARE United Arab Emirates United Arab Emirates Asia
5 ARG Argentina Argentina South America
6 ARM Armenia Armenia Asia
7 ATA Antarctica Antarctica Antarctica
8 ATF Fr. S. Antarctic Lands France Seven seas (open ocean)
9 AUS Australia Australia Oceania
10 AUT Austria Austria Europe
area pop_est pop_est_dens economy
1 652860.000 [km^2] 28400000 4.350090e+01 7. Least developed region
2 1246700.000 [km^2] 12799293 1.026654e+01 7. Least developed region
3 27400.000 [km^2] 3639453 1.328268e+02 6. Developing region
4 71252.172 [km^2] 4798491 6.734519e+01 6. Developing region
5 2736690.000 [km^2] 40913584 1.495003e+01 5. Emerging region: G20
6 28470.000 [km^2] 2967004 1.042151e+02 6. Developing region
7 12259213.973 [km^2] 3802 3.101341e-04 6. Developing region
8 7257.455 [km^2] 140 1.929051e-02 6. Developing region
9 7682300.000 [km^2] 21262641 2.767744e+00 2. Developed region: nonG7
10 82523.000 [km^2] 8210281 9.949082e+01 2. Developed region: nonG7
income_grp geometry
1 5. Low income MULTIPOLYGON (((61.21082 35...
2 3. Upper middle income MULTIPOLYGON (((16.32653 -5...
3 4. Lower middle income MULTIPOLYGON (((20.59025 41...
4 2. High income: nonOECD MULTIPOLYGON (((51.57952 24...
5 3. Upper middle income MULTIPOLYGON (((-65.5 -55.2...
6 4. Lower middle income MULTIPOLYGON (((43.58275 41...
7 2. High income: nonOECD MULTIPOLYGON (((-59.57209 -...
8 2. High income: nonOECD MULTIPOLYGON (((68.935 -48....
9 1. High income: OECD MULTIPOLYGON (((145.398 -40...
10 1. High income: OECD MULTIPOLYGON (((16.97967 48...
wmpb_World_BMI <- wmpb_World %>% left_join(big_mac,
by = c("name" = "country"))library(dplyr)
wmpb_World_BMI <- wmpb_World_BMI %>%
#filter(year == 2020) %>%
group_by(name) %>%
filter(
any(!is.na(currency_code)) &
any(!is.na(bmi_localprice)) &
any(!is.na(bmi_usd_price)) &
any(!is.na(export_usd)) &
any(!is.na(import_usd)) &
any(!is.na(GDP)) &
any(!is.na(gdp_per_capita)) &
any(!is.na(gdp_per_employed)) &
any(!is.na(inflation)) &
any(!is.na(year))
) %>%
ungroup()tmap_mode("view")
tm <- tm_shape(wmpb_World_BMI) +
tm_polygons("bmi_usd_price", id = "name", popup.vars = c("Big Mac Index(USD)" = "bmi_usd_price",
"Inflation" = "inflation",
"Population Density" = "pop_est_dens",
"GDP per Employed" = "gdp_per_employed",
"BMI local Price" = "bmi_localprice",
"Econonic Status" = "economy",
"Income" = "income_grp"),
popup.format = list(pop_est_dens = list(digits = 1), gdp_per_employed = list(digits = 0), inflation = list(digits = 1)
)) +
tm_layout()
tm#qtm(wmpb_World_BMI_2020,
#fill = "usd_price")tm_shape(wmpb_World_BMI) +
tm_fill("bmi_usd_price",
style = "quantile",
palette = "Blues",
thres.poly = 0) +
tm_facets(by="income_grp",
free.coords=TRUE,
drop.shapes=FALSE) +
tm_layout(legend.show = FALSE,
title.position = c("center", "right"),
title.size = 20) +
tm_borders(alpha = 0.5)tm_shape(wmpb_World_BMI) +
tm_fill("bmi_usd_price",
style = "quantile",
palette = "Blues",
legend.hist = TRUE,
legend.is.portrait = TRUE,
legend.hist.z = 0.1,
thres.poly = 0) +
tm_facets(by="economy",
free.coords=TRUE,
drop.shapes=FALSE) +
tm_layout(legend.show = FALSE,
title.position = c("center", "right"),
title.size = 20) +
tm_borders(alpha = 0.5) tm_layout(legend.position = c("bottom", "bottom"))$tm_layout
$tm_layout$legend.position
[1] "bottom" "bottom"
$tm_layout$style
[1] NA
attr(,"class")
[1] "tm"
tm_shape(wmpb_World_BMI) +
tm_bubbles(col = "purple",scale = 2,
size = "bmi_usd_price",
border.col = "black",
border.lwd = 3, alpha = 0.5)var <- wmpb_World_BMI$bmi_usd_price
percent <- c(0,.01,.1,.5,.9,.99,1)
quantiles <- quantile(var,percent)
print(quantiles) 0% 1% 10% 50% 90% 99% 100%
0.798722 1.251147 1.755654 3.033469 4.946211 7.064186 8.063016
get.var <- function(var_name, wmpb_World_BMI) {
# Assuming 'var_name' is a string representing the variable name to extract
v <- wmpb_World_BMI[[var_name]] %>%
as.numeric()
return(v)
}library(tmap)
tmap_mode("plot")
tmap_options(legend.width = 0.25) # Adjusts the width of the legend area
# Define the function
percentmap <- function(var_name, wmpb_World_BMI, legtitle=NA, mtitle="Percentile Map") {
percent <- c(0, .01, .1, .5, .9, .99, 1)
var <- get.var(var_name, wmpb_World_BMI)
bperc <- quantile(var, percent)
tm_shape(wmpb_World_BMI) +
tm_fill(var_name,
title = legtitle,
breaks = bperc,
palette = "Greens",
legend.hist = TRUE) +
tm_borders() +
tm_layout(main.title = mtitle,
main.title.position = "center",
legend.position = c("left", "bottom")) -> tm
print(tm)
}
# Execute the function
percentmap("bmi_usd_price", wmpb_World_BMI, legtitle = "Percentile", mtitle = "BMI USD Percentile Map")
ggplot(data = wmpb_World_BMI,
aes(x = "",
y = bmi_usd_price)) +
geom_boxplot()
geom_dotplot()geom_dotplot: binaxis = x, stackdir = up, stackratio = 1, dotsize = 1, stackgroups = FALSE, na.rm = FALSE
stat_bindot: binaxis = x, binwidth = NULL, binpositions = bygroup, method = dotdensity, origin = NULL, right = TRUE, width = 0.9, drop = FALSE, na.rm = FALSE
position_identity
boxbreaks <- function(v,mult=1.5) {
qv <- unname(quantile(v))
iqr <- qv[3] - qv[1]
upfence <- qv[3] + mult * iqr
lofence <- qv[1] - mult * iqr
# initialize break points vector
bb <- vector(mode="numeric",length=6)
# logic for lower and upper fences
if (lofence < qv[1]) { # no lower outliers
bb[1] <- lofence
bb[2] <- floor(qv[1])
} else {
bb[2] <- lofence
bb[1] <- qv[1]
}
if (upfence > qv[4]) { # no upper outliers
bb[6] <- upfence
bb[5] <- ceiling(qv[4])
} else {
bb[5] <- upfence
bb[6] <- qv[4]
}
bb <- sort(bb)
bb[2:4] <- qv[1:3]
return(bb)
}get.var <- function(vname,df) {
v <- df[vname] %>% st_set_geometry(NULL)
v <- unname(v[,1])
return(v)
}get.var <- function(vname, df) {
# This assumes 'vname' is the name of the variable to extract,
# and 'df' is the data frame or sf object without geometry.
v <- df[[vname]] %>% as.numeric() # Extract as vector and ensure it's numeric
return(v)
}
library(sf) # Assuming you're working with sf objects
library(dplyr) # For the pipe operator
var <- na.omit(get.var("bmi_usd_price", wmpb_World_BMI))
bb <- boxbreaks(var)
print(bb)[1] -2.553398 0.798722 2.324783 3.033469 5.000000 6.385589
boxmap <- function(vnam, df,
legtitle=NA,
mtitle="Box Map",
mult=1.5){
var <- get.var(vnam,df)
bb <- boxbreaks(var)
tm_shape(df) +
tm_polygons() +
tm_shape(df) +
tm_fill(vnam,title=legtitle,
breaks=bb,
palette="Oranges",
labels = c("lower outlier",
"< 25%",
"25% - 50%",
"50% - 75%",
"> 75%",
"upper outlier")) +
tm_borders() +
tm_layout(main.title = mtitle,
title.position = c("left",
"top"))
}tmap_mode("view")
boxmap("bmi_usd_price", wmpb_World_BMI)+
tm_facets(by="continent",
free.coords=TRUE,
drop.shapes=FALSE)